#include "AnaDef.h"

const Double_t EPSILON = 1e-12;

TGraphAsymmErrors * getRatio(TGraphAsymmErrors * grUp, TGraphAsymmErrors * grdown)
{
  TGraphAsymmErrors * grr = new TGraphAsymmErrors;

  for(Int_t ii = 0; ii< grUp->GetN(); ii++){
    const Double_t xx = grUp->GetX()[ii];
    if(fabs(xx-grdown->GetX()[ii])>EPSILON){
      printf("bad x %d %e %e\n", ii, xx, grdown->GetX()[ii]); exit(1);
    }

    const Double_t exl = grUp->GetEXlow()[ii];
    const Double_t exh = grUp->GetEXhigh()[ii];

    const Double_t yUp = grUp->GetY()[ii];
    const Double_t elUp = grUp->GetEYlow()[ii];
    const Double_t ehUp = grUp->GetEYhigh()[ii];

    const Double_t ydown = grdown->GetY()[ii];
    const Double_t eldown = grdown->GetEYlow()[ii];
    const Double_t ehdown = grdown->GetEYhigh()[ii];

    const Double_t yy = yUp/(ydown+EPSILON);
    const Double_t el = yy*TMath::Sqrt(TMath::Power(elUp/yUp,2)+TMath::Power(eldown/ydown,2));
    const Double_t eh = yy*TMath::Sqrt(TMath::Power(ehUp/yUp,2)+TMath::Power(ehdown/ydown,2));

    grr->SetPoint(ii, xx, yy);
    grr->SetPointError(ii, exl, exh, el, eh);
  }

  return grr;
}

void yieldRatio(const TString fin)
{
  TFile * ff= new TFile(fin);
  if(!ff){
    printf("no ff! %s\n", fin.Data()); exit(1);
  }

  TGraphAsymmErrors * grstat[]={0x0,0x0,0x0,0x0};
  TGraphAsymmErrors * grsys[]={0x0,0x0,0x0,0x0};
  for(Int_t ii=0; ii<4; ii++){
    grstat[ii]=(TGraphAsymmErrors*) ff->Get(Form("grstat%d", ii));
    grsys[ii] =(TGraphAsymmErrors*) ff->Get(Form("grsys%d", ii));
    if(!grstat[ii] || !grsys[ii]){
      printf("grstat or sys bad %d\n", ii);ff->ls(); exit(1);
    }
  }

  TGraphAsymmErrors * outStat[]={0x0,0x0,0x0,0x0};
  TGraphAsymmErrors * outSys[]={0x0,0x0,0x0,0x0};
  for(Int_t ii=0; ii<4; ii++){
    outStat[ii] = getRatio(grstat[ii],grstat[1]);
    outStat[ii]->SetName(Form("grstat%d", ii));

    outSys[ii] = getRatio(grsys[ii],grsys[1]);
    outSys[ii]->SetName(Form("grsys%d", ii));
  }

  TString fname(fin);
  while(fname.Contains("/")){
    fname=fname(fname.First('/')+1,fname.Length());
  }
  
  fname.Prepend("fromSpe_");

  TFile *fout =new TFile(fname,"recreate");
  for(Int_t ii=0; ii<4; ii++){
    outStat[ii]->Write();
    outSys[ii]->Write();
  }
  fout->Save();
  fout->Close();
  delete fout;
}


int main(int argc, char* argv[])
{
  for(Int_t  ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  if(argc!=2){
    printf("argc!=2\n"); return 1;
  }

  yieldRatio(argv[1]);

  return 0;
}
